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ABSTRACT 

A numerical study of the long-time evolution of incompressible Navier-Stokes turbulence forced 
at a single long- wavelength Fourier mode, i.e., a Kolmogorov flow, has been completed. The 
boundary conditions are periodic in three dimensions and the forcing is effected by imposing 
a steady, two-dimensional, sinusoidal shear velocity which is directed along the x- direction and 
varies along the ^-direction. A comparison with experimental data shows agreement with mea- 
sured cross-correlations of the turbulent velocity components which lie in the mean-flow plane. 
A statistical analysis reveals that the shear-driven turbulence studied here has significant spec- 
tral anisotropy which increases with wave number. 


a This research was supported by the National Aeronautics and Space Administraton under NASA Contract 
No. NAS 1-19480 while the second author was in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23681-0001. 
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I. Introduction 


The presence of an unstable mean shear (Le., a nonzero time- averaged velocity gradient) is the essential 
source of turbulent fluctuations in any given fluid dynamical system. The mean shear is maintained by 
suitable external forces and, if the system contains a region of fully developed turbulent flow which exists 
long enough for study, then meaningful statistical statements can be made concerning it. For example, the 
characteristics of statistically stationary flow through pipes and channels [1, 2, 3, 4] and over flat plates [5, 6] 
have been carefully observed experimentally for over a century. 

While the laminar pipe flow profile is parabolic, an even simpler laminar profile is a linear one. Thus, in 
the last several decades there have been a number of physical experiments measuring the characteristics of 
incompressible fluid turbulence in the presence of a steady linear velocity gradient [7, 8, 9, 10] (with regard 
to the results of [8], however, see [11]). This basic flow has also been simulated numerically, although the 
numerical procedure is rather complicated, since difficulties arise in treating the turbulent flow as homoge- 
neous, (Le., with periodic boundary conditions) while simultaneously imposing a non-periodic linear shear 
profile to drive the flow [12, 13, 14, 15]. 

Here, another basic flow is utilized to examine forced, dissipative turbulence: An erstwhile homogeneous 
flow forced by a sinusoidal shear profile. This steady sinusoidal shear has a wavelength equal to an edge-length 
of the periodic box, where the imposed velocity is always in the ^-direction but varies along the 2 -direction. 
The resulting motion, which can be considered as homogeneous in each a^-plane, but not in the 2 -direction, 
is a Kolmogorov flow. In previous work, transition and turbulence of two-dimensional Kolmogorov flows have 
been discussed [16, 17, 18, 19, 20, 21, 22] and a study of three-dimensional Kolmogorov flow in the presence of 
hyperviscosity was reported [23]; here, three-dimensional Kolmogorov flow turbulence with normal viscosity 
is numerically simulated and examined. 

Near the zeroes of the sinusoidal shear velocity the profile is approximately linear, while near the extrema 
of the imposed velocity the profile is approximately parabolic. Heuristically, the basic Kolmogorov flow can 
be thought of as a combination of a linear homogeneous shear profile and a laminar parabolic channel 
profile. As with these profiles ( i.e Couette and Poiseuille flow [24]), an imposed sinusoidal shear velocity 
of sufficiently large amplitude is also unstable to perturbations [19]. Here, the amplitude is large enough so 
that instability occurs and turbulence results. 


1 



In brief, the purpose of this paper is to report on Fourier-method simulations of three-dimensional 
Kolmogorov turbulence. This will provide a numerical flow history with which comparisons to physical 
experiments can be made. This numerical record will then allow some aspects of physical and spectral 
anisotropy to be illustrated. 

In the next section, the basic equations and qualitative features of their solution are presented. Fol- 
lowing this, the numerical procedure is briefly discussed and numerical results are compared with pertinent 
experimental ones. Finally, a discussion and conclusions are given. 

II. Basic Equations 

In vorticity form, the equations which describe incompressible Navier-Stokes flow are 

d t fi = V x (V x f?) + i/V 2 ft + V x f (1) 

Here, the total fluid velocity is V, the vorticity is ft — V x V, and f is the external forcing applied to the 
fluid; for incompressibility, we also have V * V = 0. Equation (1) is non-dimensional, so that v = Rj 1 , 
where R e is the Reynolds number; additionally, the density has been set to unity. If the fluid is driven by 
an imposed velocity field V 0 (where V * V 0 = 0, with associated vorticity ft 0 = V x V 0 ), the total velocity 
and the total associated vorticity can be written as 

. V = V 0 -j- u and ft = ft 0 4- u? (2) 

where u is the velocity (satisfying V • u = 0) and w — V x u is the vorticity which arise from perturbations 
to the Kolmogorov flow: 

V 0 .== xA sin z, and ft 0 = yA cos z . (3) 

(When needed, the unit vectors x, y, and z, in the z-, y-, and z-directions, respectively, will be used). Placing 
(2) and (3) into (1) yields the basic equation of this paper: 

dt& ~ V x [(u + V 0 ) x (a? + f? 0 )] T z^V 2 cj. (4) 

In order to produce (4), the external force in (1) is f = x vA sin z; also, using (3), note that V x (V 0 x ft 0 ) — 0. 
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III. Qualitative Analysis 


The model physical space in which the forced, dissipative turbulent flow is to occur will be a three-dimensional 
cube on which periodic boundary conditions are imposed. A qualitative picture of solutions to (4) can then 
be gained by using a discrete Fourier transformation to approximate these partial differential equations by 
a much larger, but finite, set of coupled ordinary differential equations. To do this, we use 

u(x,£) = u(k,t)e^ x and w(x,t) = ^ a?(k, ( 5 ) 

|k|<iv /2 |k«iv /2 

Here, i — and the sum (and each similar sum appearing henceforth) is over all k such that k = jkj < 

| AT, i.e., the Fourier series is ‘isotropically truncated 9 to ensure no inherent spectral bias; N is the number 
of points in each of the three spatial dimensions. The flow field is thus represented by commensurate sets 
of values on discrete points in physical space and in Fourier space: x — (xi,yj,Zk) ~ 7 f (h j i k) ) i,j y k G 
{0, 1, • ■ TV— 1}, and k ■= (/, m, n), /, m, n G 1, * * *, — 1,0, 1, * • •, |iV} with k = y/l 2 -f m 2 -f n 2 < |i\T. 

The velocity and vorticity coefficients are related by u(k) = ik ~ 2 k x u>(k) and o;(k) = ik x u(k). Since 
the physical space velocity and vorticity fields u(x) and u?(x) are real, their Fourier coefficients satisfy 
u(k) ■== u*(— k) and u?(k) = u>*(— k), where signifies complex conjugation (and explicit time-dependence 
in the argument has been dropped for brevity). 

Putting (4) and (5) together and using w(k) = du(k)/dt leads to 

d?(k) ,== «k x ^2 [u(p) -f V 0 (p)] x [o?(k - p) + SI 0 { k - p)] - i/k 2 co( k). ( 6 ) 

IPKJV/2 

Here the Fourier coefficients associated with (3) are 

V 0 (±k) = ±t^8(z qp k)x and f2 0 (dbk) = -y 8(z 3 = k)y (7) 

where <5(p) = 0 if p ^ 0 and <5(p) = 1 if p = 0 [i.e., 6 (p) is the Kronecker delta function]. 

A complete qualitative analysis requires determining the critical -points of ( 6 ), i.e., the zeros of the right- 
hand-side. Although a full discussion is beyond the scope of this work, one critical point of ( 6 ) occurs when 
all u;(k) = 0 , that is, at the origin of the phase space of the independent real and imaginary parts of the 
Fourier coefficients. If we linearize ( 6 ) about this point, we have 

w(k) = R + (k) * w(k - z) + R-(k) • w(k + z) - vk 2 u>( k). ( 8 ) 


3 



The tensors R ± (k) are real quantities; their components are 



A , 

^ e a7P fc 7 


(k y S p( 3 ± Spzfyy) 


|kTz 


$12 


± S pz$Py) 


( 9 ) 


Here, repeated indices are summed over and € ayp is the completely antisymmetric tensor with t xyz = 1, 
The coupled set of equations (8) presents a complicated eigenvalue problem which must be addressed to 
fully understand the nature of the critical point at the origin of the phase space. This will not be done here; 
instead, consider the following two-dimensional system, which will illustrate some general features of (8): 


p = Bq — vp 

q = Cp—vq (10) 

The system matrix implicit on the right-hand-side has eigenvalues X± = —v ± y/'BC , so that the two 

independent solutions to (10) have time dependence ~ exp(A ±t). Assuming BC ps ±A 2 with A > 0, there 

are three cases: First, let BC = —A 2 ; then the two-dimensional system point will merely spiral into the 
origin, so that p, q — ► 0 as t — ► oo. Second, let BC = A 2 with 0 < A < v\ in this case, the system point will 
head directly to the origin without spiraling (unless A := v, in which event it moves to a point offset from 
the origin). In the third case, let BC ~ A 2 with v < A; here, in general, the system point will diverge away 
from the origin exponentially with time. 

Relating this simplified system behavior to that of (8), it is clear that for A sufficiently small in (9), the 
phase point associated with (8) will return to the origin for any small perturbation so that the origin is a 
stable fixed point - an attractor. However, as A is increased, there will come a time at which the phase 
point diverges in at least one direction in phase space; the associated mode will continue to grow until the 
system behavior is no longer linear; energy now leaves the growing mode not by dissipation as with the other 
modes, but by nonlinear transfer to these other modes, where it is eventually dissipated. Increasing A even 
further leads to more modes acting as energy sources (and fewer as sinks). The neighborhood of the critical 
point at the origin thus gains a multidimensional saddle point structure. (This is essentially the transition 
to turbulence described by Landau & Lifshitz [25], and for two-dimensional Kolmogorov flow, it has been 
demonstrated by She [19]). 

Suppose that A is sufficiently large so that u(k) = 0 is an unstable equilibrium; then, small perturbations 
will grow exponentially for some u(k). As the flow develops, the u(k) will generally behave as a collection 
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of random variables, i.e., u(k) will have two parts, a time-averaged part u(k) and a fluctuating part v(k): 

u(k) = u(k) + v(k) u(k) = ~ j u(k) dt (11) 

where k * u(k) := k * v(k) = 0 (and the vorticity associated with v(k) will be denoted w(k) = ik x v(k)). In 
this case, u(k) may be another critial point. (Note that in [23] the mean flow is defined as the instantaneous 
flow associated with the mode u(k^), where k / = (0, 0, 1) here; this obviates a direct comparison of results). 
Alternatively, u(k) may not have a simple attractor (a fixed point), but may have a limit cycle (a closed 
path in phase space), or even a strange attractor (an open path which is confined to a subspace of phase 
space which is often of fractal dimension [26]). 

One way to determine which possibility actually occurs is to proceed numerically. This has, in fact, already 
been done for two-dimensional Kolmogorov flow transition and turbulence [19] and for three-dimensional, 
hyperviscous Kolmogorov flow [23]. For sufficiently large Reynolds number (R e = v~ x > 30 in [19] and 
R e > 13 in [23]), turbulence occurs; here v~ x = 2000 (although the microscale Reynolds number R\ « 30). 
Thus, only the turbulent flow regime is examined here (although the nature of the attractor will not be 
considered in any further detail). 

IV Results 

The basic equations in Fourier space (6) were solved numerically using a pseudospectral method, with N 
in (5) equal to 64. To integrate (6) forward in time, the dissipative term was treated implicitly using a 
backward Euler scheme, while the nonlinear terms were handled using a 'third-order partially corrected 
Adams-Bashforth 5 method [27]. The Fourier coefficients of the nonlinear terms were determined by first 
transforming the various factors into physical space, then forming the point-by-point products, and finally 
transforming these products back into Fourier space [28]. 

The time-step was At = .0025. The initial turbulence spectrum was |u(k)| ~ & 2 exp—& 2 /4, with the 
u(k) having random phase; the total initial turbulent energy was 0.0025. There was an initial transient 
period during which the value of A in (7) was increased to A = 1.5 (from A — 0.1), and the viscosity in (6) 
was decreased to v = 0.0005 (from v — 0.001). After a quasi-steady state appeared to become established, 
the values of all the Fourier coefficients u(k) were saved for all integer values of t from t = 150 to t = 650, 
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inclusive, for a total of 501 time slices upon which to perform statistical evaluation. 

During the data collection time, t = 150 to t = 650, the parameter values were fixed at A = 1.5 and 
v = 0.0005. Denoting the volume average of any quantity Q by < Q >, the root-mean-square (rms) shear 
value can be defined as S =< ft 2 > 1 ^ 2 ; using (3) and A ■= 1.5, we have S & 1. Thus the time variable t used 
here is equivalent to the scaled time variable used in [15]. 

In Figures 1 and 2, some of the global quantities for the simulation between t = 150 and t == 650 are 
shown. In Figure 1, the energy and the enstrophy are presented: the energy is defined as E = \ < v 2 > 
and the enstrophy as Q = | < w 2 >. Thus, the rms velocity fluctuation appears to be around 0.08 
and the rms vorticity fluctuation appears to be around 0.22. In Figure 2, the following quantities are 
given: the average or rms wave number, k avg — ^/Q/E = (/ k 2 E(k)dk/ f E(k)dk) 1 ^ 2 ; the dissipation wave 
number, k& = (2 Q) 1 / 4 //" 1 / 2 ; the microscale Reynolds number, R\ == y/2E/{vk a vg) ~ ( kz)/k aV g ) 2 . Since 
kjg < N/2 = 32, the simulation is apparently well resolved. This is also borne out by the average energy 
spectrum shown in Figure 3, which is the average of the 501 individual energy spectra for each of the time 
slices between t == 150 and t = 650. The average energy spectrum is approximately exponential, especially 
in the smaller scales, with E(k) ~ e“°‘ 26fc . 

Mean Kolmogorov flow is strongly anisotropic since it is driven by the imposed velocity (3), The associated 
turbulence may also be anisotropic; a detailed examination can be conducted by considering the spectral 
anisotropy of the turbulence. Two measures of this are the ‘Kolmogorov anisotropy’ and another we will 
term the ‘& 2 -anisotropy\ The Kolmogorov anisotropy K a p(k) (a, /? € {x, y, z}) is defined as follows: 


KafiW = 


Sgp(k) 

S XX (k) + S y y(k) + S ZZ (kY 


S*p(k) = ^(k)^(k) 

|k|=i 


( 12 ) 


while the A; 2 -anisotropy C a (k) (a £ {x,y, z}) is defined as: 


C a (k) = 


Dajk) 

Dx(k) + D y (k) + D z {k) ’ 


D a (k)= Yj ^l v (k)| 2 - 
|k|=fc 


(13) 


The Kolmogorov anisotropy measures the spectral distribution of turbulent energy amongst the different 
components of the velocity, and thus quantifies the behavior of spectral transfer between the velocity compo- 
nents as a function of increasing wave number. On the other hand, & 2 -anisotropy measures the geometrical 
shape of the total spectral energy distribution, z.e., to what extent it has spherical, spheroidal, or ellipsoidal 
symmetry. 
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In Figures 4, 5, and 6, the results of time-averaging from t = 150 to i = 650 are presented: In Figure 4, 
the diagonal components of K & p(k) and in Figure 5, the off-diagonal components of K a p{k) are presented; 
in Figure 6, the components of C a {k) are given. It is clear that spectral anisotropy persists and, in fact, 
increases as k increases, for the diagonal components of I< a p(k) and for the components of C a (k). (Isotropy 
occurs when these components are all equal to |.) 

As already mentioned, a number of experimental results have appeared for the case of turbulence in the 
presence of a linear shear (V^ = zx) [7, 8, 9, 10]. The sinusoidal shear (3) used here approximates a linear 
shear around the values z = i r and z = 2w = 0 mod 27 r, so that there may be a basis for comparison. To 
effect a comparison, the following xy-plane-averaged cross-correlations Q a p and stress anisotropy vector B a , 
related to the Reynolds stress tensor, were averaged and symmetrized in physical space: 

Qccp(z) = / < Va - P = ==> a P € {*!/> v z > zx }’> B °( z ) = < ^ — -l, a € {x, y, z}. (14) 

*J< v% > z < vj > z E * 3 

The brackets < * • * >* indicate an average over all points in the £?/- plane in physical space corresponding to 

each of the N values of z; E z = | < v 2 > z is the kinetic energy in a given z-plane. These planar averages 

were taken for each of the 501 time slices between t = 150 and t = 650, and then averaged together. These 

time-averages were further averaged by symmetrizing : 

Qap{z) = jiQap^ + Qapi^TT- z)-Q a p(w- z)-Qai}(ir+z)] 

B a (z) = + B a ( 2ir - z) + B a (it - z) + B a { n + z)j (15) 

where it must be remembered that all flow variable are periodic in z (so that z z mod 27r). 

The results of this symmetrization are shown in Figure 7 for Q zx (z ) and in Figures 8, 9, and 10 for 
B a (z). (Components Q xy (z) and Q yz {z) were essentially zero for all z, and will not be shown.) The value 
\Qzx{z)\ & 0.45 at z = 0 and tt correlates extremely well with the existing experimental results, which are 
summarized in Table 1 (which itself draws from another summary appearing in [33]). Although the B a (z) 
values do not match up very well, there is also a substantial variation amongst the various experimental 
values, as Table 1 again shows. A possible explanation for large difference between values of the B a {z) 
for the flow examined here (sinusoidal shear) and those derived from linear shear experiments may lie in 
the enhanced mixing of the sinusoidal shear, which tends to equalize the average magnitudes of v x and v z . 
Additionally, the close match of Q zx may indicate a more universal measure of shear-driven turbulence than 
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the Bq,(z) can provide. 

V Discussion 

The results shown in Figures 4 and 6 do not appear consistent with the assumption of isotropy on which the 
Kolmogorov theory of turbulence, which predicts that the energy spectrum in the ‘inertial range 5 (kj « 
k « fep) has the form E(k ) ~ fc” 5/3 , is partly based [29]. However, another assumption in the Kolmogorov 
theory, that of homogeneity, is also not met here; the presence of sinusoidal forcing means that homogeneity 
can only be assumed to exist in two dimensions. Additionally, the relatively low Reynolds number and 
consequent lack of an inertial range here (and in any practical direct numerical simulation) indicates the 
absence of the very spectral regime for which homogeneity is a useful concept. 

Nevertheless, it has also been argued [30] “that in uniform shear flow near equilibrium, local [physical] 
isotropy can never constitute a systematic approximation, even in the limit of infinite Reynolds number.” 
This view is, in fact, borne out by physical experiments over a wide variety of Reynolds numbers [2, 5, 31], 
where the relative strengths of fluctuations in different coordinate direction are consistently anisotropic (as 
Table 1 shows). In spite of the anisotropy, the &“ 5 / 3 spectrum is still seen [2], particularly in the streamwise 
fluctuations. 

The Kolmogorov spectrum E(k) ~ fc~ 5 / 3 , which is well established experimentally in many types of 
flows [32], may depend more on its dimensional consistency and the assumption of power-law behavior, 
than on other prior assumptions (such as isotropy at the small scales, i.e., at large values of wave number). 
Nevertheless, Figures 4 and 6 unequivocally show anisotropic spectral behavior with increasing wave number, 
while Figures 8,9, and 10 slow only slight anisotropy in physical space due to the bias of these physical-space 
measures of anisotropy toward the relatively isotropic large scales. Whether or not the spectral anisotropy 
is maintained with increasing numerical resolution is an open question. 

VI Conclusions 

In this paper, the results of a numerical simulation of three-dimensional Kolmogorov flow were presented. 
This simulation was done on a 64 3 grid, and run for a relatively long period: T tot — 500 (2 x 10 5 time steps 
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at At 5= .0025), after a period of initial transience had passed (as already mentioned, this time is equivalent 
to the rms-shear-scaled time of [15]). The forcing wave number was the lowest possible ( kj •= 1), in contrast 
to that of some previous two-dimensional studies [19, 22], where kj = 8 was used. The numerical simulation 
produced a flow history consisting of 501 sets of complete Fourier coefficients, each set effectively a time-slice 
of the dynamic evolution, for every whole value of simulation time between t = 150 and t = 650. In analyzing 
the flow history, it was evident that the flow was in a near-equilibrium state (although inherent fluctuations 
will always be present). 

One interesting result of the simulation presented here concerned the (normalized) Reynolds stress tensor: 
The ‘cross-correlation’ Q zx ss 0.45, in keeping with experimental linear shear results and providing some 
validation of the numerical simulation. However, the stress anisotropy vector components were more nearly 
equal than experimental linear shear results portrayed, though it is clear that the B a from different sources 
had a large amount of variation, particularly in comparison to the Q zx (as Table I shows). This corroborates, 
to some extent, those turbulence modelling assumptions which presume that the B a are essentially equal in 
fully turbulent flows.. 

The second interesting result of the simulation was that the turbulent fluctuations in the long-wavelength 
Kolmogorov flow appeared to have significant spectral anisotropy which increased with wave number. This 
anisotropy presumably resulted from the presence of forcing, as unforced turbulence was seen to be essen- 
tially isotropic at all wave numbers [34]. The importance of this result is that the Kolmogorov flow studied 
here is analogous to realistic shear flows, such as those through channels and over bounding surfaces (z.e., in 
boundary layers). If, in such cases, the flow also shows spectral anisotropy (and increasingly so) at higher 
and higher wave numbers, then this basic finding represents a new factor for incorporation into turbulence 
models. However, the assumption of the isotropy of those physical quantities which are primarily determined 
by the large scales of the turbulent motion is not invalidated as a useful approximation by the present results. 

Acknowledgements . We would- like to thank G. Erlebacher, S. Girimaji, and Y. Zhou of ICASE and B. 
Younis of City University, London, for their useful comments. 
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Table 1. Comparison with Experimental Results. 


Source 

~~Qzx 

B x 

By 

B z 

Linear shear ([33], Table 1) 

0.48 

0.30 

-0.12 

-0.18 

Near-wall turb. ([33], Table 1) 

0.45 

0.51 

-0.09 

-0.42 

Linear shear ([9], Table 4) 

0.45 

0.39 

-0.10 

-0.29 

Sinusoidal shear (this paper) 

0.45 

0.08 

-0.07 

-0.01 
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Figure 1. Turbulent energy E and enstrophy 0 versus time. 
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Figure 2. Average wave number k avg , dissipation wave number &d, and microscale Reynolds 
number R\ versus time. 
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Figure 3. Averaged turbulent energy spectrum -E(fc) versus wave number magnitude k. 
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Figure 4. Averaged diagonal Kolmogorov anisotropy components versus wave number mag- 
nitude. 
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Figure 7. Averaged cross-correlation Q zx (z) versus z. 









Figure 10. Averaged stress anisotropy component B z (z) 


versus 2 . 
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